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Abstract 

Patterns of different symmetries may arise after solution to reaction- 
diffusion equations. Hexagonal arrays, layers and their perturbations 
are observed in different models after numerical solution to the cor¬ 
responding initial-boundary value problems. We demonstrate an in¬ 
timate connection between pattern formations and optimal random 
packing on the plane. The main study is based on the following two 
points. First, the diffusive flux in reaction-diffusion systems is ap¬ 
proximated by piecewise linear functions in the framework of struc¬ 
tural approximations. This leads to a discrete network approximation 
of the considered continuous problem. Second, the discrete energy 
minimization yields optimal random packing of the domains (disks) 
in the representative cell. Therefore, the general problem of pattern 
formations based on the reaction-diffusion equations is reduced to the 
geometric problem of random packing. It is demonstrated that all 
random packings can be divided onto classes associated with classes 
of isomorphic graphs obtained form the Delaunay triangulation. The 
unique optimal solution is constructed in each class of the random 
packings. If the number of disks per representative cell is finite, the 
number of classes of isomorphic graphs, hence, the number of optimal 
packings is also finite. 


1 Introduction 

The Turing mechanism for reaction-diffusion equations models biological and 
chemical pattern formations. This approach was widely discussed in litera¬ 
ture and supported by many numerical examples (see the recent books 
[3] and many works cited therein). Patterns of different symmetries may arise 
after solution to reaction-diffusion equations. Hexagonal arrays, layers and 
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their perturbations are observed in different models after numerical solution 
to the corresponding initial-boundary value problems for nonlinear partial 
differential equations. However, these models do not answer the question, 
why the most frequently observed patterns are close to the optimal packing 
structures. Why do the hexagonal array arise? One can see, for instance, 
that a resulting structure can be the hexagonal array disturbed by pentagon 
inclusions. Is it related to a model approximation or to an inherent feature 
of pattern formations? 

In the present paper, we try to answer the above questions to demon¬ 
strate an intimate connection between pattern formations and optimal ran¬ 
dom packing on the plane. The main study is based on the following two 
points. First, the diffusive flux in react ion-diffusion systems is approximated 
by piecewise linear functions in the framework of structural approximations 
12 ], This leads to a discrete network approximation of the considered 
continuous problem. Second, the discrete energy minimization yields opti¬ 
mal random packing of the domains in the representative cell. The packed 
domains are approximated by equal disks. These approach is described in 
the bulk of the paper. 

Packing problems refers to geometrical optimization problems [10] . In the 
present paper, we consider the optimal packing of disks on the plane in the 
random statement htted to the description of pattern formations. Optimal 
packing in the classic deterministic statement is attained for the hexagonal 
array when the packing concentration holds uni. Computer simulations 
demonstrate that random packing have a lower density and depends on the 
protocol of the random packing. 

It is shown in SecJS] that pattern formations lead to the optimal ran¬ 
dom packing problem. In the present paper, this problem is resolved by 
introduction of the equivalence classes of graphs obtained form the Delaunay 
triangulation associated to packings. The justihcation of such an approach is 
based on the observation that solution to the physical problem of the optimal 
diffusion implies solution to the geometrical problem of the packing disks [7|. 
The unique optimal solution is constructed in each class of the random pack¬ 
ings. If the number of disks per representative cell is finite, the number of 
classes of isomorphic graphs, hence, the number of optimal packings is also 
hnite. 

The proposed method to study pattern formations is based on the inves¬ 
tigation of graph structures by analytical and numerical methods without 
treatment of PDF. 
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2 Structural approximation 

The Turing mechanism can create temporally stable and spatially non-homogeneous 
structures. In order to present the main idea of the structural approximation 
we consider ID Schnakenberg system [1[ p. 156]. A typical dependence of 
the inhibitor on the spatial variable is displayed in FigJT^. It is assumed that 
such a dependence can be approximated by a piecewise linear function as 
shown in FigOj. The solution of the continuous reaction-diffusion equations 



a) 



Figure 1: a) Dependence of the inhibitor on the spatial variable, b) Piecewise 
linear approximation of the inhibitor on a smaller interval (dashed line). The 
maxima are approximated by segments and Dj (disks in 2D) and the 
minima by points (segments in 2D). 

is approximated by the discrete diffusion model with the constant diffusion 
fluxes (derivatives of the linear approximations) between the extrema of the 
potential. 

A similar approximation can be extended to multidimensional reaction- 
diffusion equations [S]. In the present paper, we deal with 2D double periodic 
structures. Let ei = (ei,0) and 62 = ( 621 , 622 ) be the translation vectors of 
the lattice Q = {/iCi -|- ^262 : Ij G Z} where Z denotes the set of integer 
numbers. Consider the periodic representative cell 

Qo = {x = tiGi ^ 2 ^ 2 }, 0 < tj < 1}. 

For simplicity, we approximate the places of maximal diffusion potential by 
equal disks Di {i = 1,2,N) of radius r centered at the set of points A = 
(ai, a 2 ,..., slm) displayed in Figj2j The maxima of the diffusion potential are 
approximated by disks and the minima by segment^. Every line segment 

^It is natural to introduce the disk approximations also for minima. But in this case 
we shall obtain two types of disks that complicates the one-type-disks model constructed 
below. Because the diffusion flux will occur between the disks of different types and vanish 
between the disks of the same type. 
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Figure 2 : 2 D approximation of the inhibitor. The diffusion potential is ap¬ 
proximated by appropriate constants in disks and the diffusion flux between 
the disks by linear functions along the edges of the Delaunay triangulation. 

Pkj is perpendicular to the segment (a^jUj), its length holds \Pkj\ = 2r and 
it is divided onto equal parts by (a^, a^). 

It is convenient to introduce new distance (metric) as follows. Two points 
a, b G are identihed if their difference a — b = liGi -|- (262 belongs to the 
lattice Q. Hence, the classic flat torus topology with the opposite sides 
welded is introduced on the cell Qq. The distance ||a — b|| between two 
points a, b G Qo is introduced as 

lla —b||:= min la — b-|-/iCi-|-/2e2| , (2.1) 

where the modulus means the Euclidean distance in between the points 

a and b. 

Construct the double periodic Voronoi diagram and the Delaunay triangu¬ 
lation corresponding to the set A on the torus Qq = h^i + ^^2)■ 

The edges of the Delaunay triangulation E correspond to linear approxima¬ 
tions of the diffusion flux between disks. The Delauney triangulation of the 
vertices A consists of straight lines connecting by pairs points of A belong¬ 
ing to neighbor Voronoi region^. Let the neighborhood relation between two 

^The terms the Delaunay triangulation and graph used in this paper are slightly dif¬ 
ferent from the commonly used notations in degenerate cases. For example, consider a 
square and its four vertices. The traditional Delaunay triangulation has four sides of the 
square and one of the diagonals. In our approach, the Delaunay graph has only four sides. 
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vertexes be denoted by ~ or shortly j ~ k. We call the constructed 
double periodic graph (A, E) by the Delaunay graph. 

Two graphs are called isomorphic if they contain the same number of 
vertices connected in the same way. One of the most important notation of 
the persent paper is the class of graphs Q = G{a,e) isomorphic to the given 
graph {A,E). 

Let u = {ui,U2, ■ ■ ■ ,U]sf) denote the vector whose components are the 
maximal diffusion potentials in the corresponding disks. The discrete network 
model for densely packed disks 13. 0. 13 is based on the fact that the 
diffusion flux is concentrated in the necks between closely spaced inclusions 
having different potentials. In our model, closely spaced inclusions means the 
chain disk-segment-disk Pkj Dj) displayed in Fig{T]D. For two 

neighbor disks Dj with centers a^, and a segment Pkj between them 
the relative interparticle flux ^^da^ — aj|) can be approximated by Keller’s 
formula [5] 



( 2 . 2 ) 


where 5kj = ||afc — aj|| — 2r denotes the gap between the neighbor disks. 
Keller’s formula fl2.2p was deduced for the linear local diffusion flux between 
the neighbor disks that is agree with our approximation. Introduce the des¬ 
ignation 



(2.3) 


where j ~ /c means that the vertices aj and a^ are connected. Following [2], 
[6] introduce the functional associated to the discrete energy 



(2.4) 


where Uk + Uj is the variation of the diffusion potential along the chain 

Dk Pkj Dj- 

3 Optimal random packing 

Consider the minimization problem 


T(u) = min i?(u, A) = min 


A A 


g{\\ak - S-j\\)iuk + Ujf. (3.1) 
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The function g{x) = as a convex function satisfies Jensen’s inequality 


M 


M 


^Pi9{xi) > g ( ^PiXi , , 

i=l 


(3.2) 


i=l 


where the sum of positive numbers Pi is equal to unity. Equality holds if and 
only if all Xj are equal. Let the sum from fl2.4l) is arranged in such a way that 
Xi = ||afc —ajil and Pi = ^{uk + Ujb, where U = Y^^^\uk + Ujb. Application 
of 03.21) to 02.41) yields 


- ^j\ 


i||)(^fc “1“ ^j) ^ Ug ^ ^ (^/l “1“ 11^/^ I 

Holder’s inequality states that for non-negative a* and bi 


(3.3) 


M 


M 


2 / M 




(3.4) 


2 = 1 


. 2=1 


. 2=1 


This implies that 

\2ll II ^ ^ II 

/ j ||afc II — / j T ^j) / ^ ll^fc 


(3.5) 


The function g{x) decreases, hence 03.3p and 03.5p give 


3{\\^^-^^\\){uk+Ujf >g[Yj + ||afc - a 


(5), 


siG) 


U 


u 


(3.6) 

The minimum of the right hand part of 03.61) on A is achieved independently 
on Uk for max A h{A) where 


^(A) = l|at - a 


(3.7) 


Lemma 3.1 ([S])- For any fixed class G(a, E), every local maximizer of h{A) 
is the global maximizer which fulfils the system of linear algebraic equations 


^ jr^k ^ £=1,2 


(3.8) 


where Ski can take the values 0,±1,±2 in accordance with the class G{a,e)- 
The system 03. 8 p has always a unique solution up to an additive arbitrary 
constant vector. 
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Equations fl3.8p describe the stationary points of the functional fl3.7p ob¬ 
tained by its differentiation on (A: = 1, 2,, N) 

- afc) = 0. (3.9) 

j^k 

Here, the congruence relation a = b means that a — b = liei + ^ 2^2 for some 
integer li^ 2 - Therefore, a point a on the torus Qq is associated to the inhnite 
set of points {a -|- liSi + I 2 G 2 , h ,2 £ on the plane M^. We now rewrite 
equation fl3.9p on the torus as an equation on the plane for a hxed point 
£ Qo- Consider a points a)- G neighboring to a^, i.e., j ~ A: in a graph 
(A, E) G 0(A,E)- The point a), is congruent to a point a^ G Qo- The graph 
(A, E) corresponds to the Voronoi tessellation, hence, a'- belongs to Qo or to 
neighbor cells Qo±ei, (5o±e2, Qo±ei±e 2 . Therefore, a'- = SLj+lijkei+l 2 jke 2 , 
where hjk and l 2 jk can be equal only to 0, ±1. Then, equations 03.91) can be 
written in the form 03.81) where 

■Slfc ^ ^ 0jky ^2k ^ ^ ^2jk- (3.10) 

j'^k j'^k 

One can see that the sum of all equations 03.8p gives an identity, hence, 
they are linearly dependent. Moreover, if A = (ai, a 2 ,..., ajv) is a solution 
of 03.8p . then (ai -|- c, a 2 -|- c,..., a^r -|- c) is also a solution of 03.8p for any 
c G M^. Let the point uat be arbitrarily hxed. Then, ai, a 2 ,..., aAr_i can be 
found from the uniquely solvable system of linear algebraic equations 

A; = 1 , 2 ,..., A - 1 . (3.11) 

^ j^k ^ 1=1,2 

It is worth noting that the system 03.lip can be decomposed onto two inde¬ 
pendent systems of scalar equations on the hrst and second coordinates of 
the points ai, a 2 ,..., aAr_i. 


4 Conclusion and numerical examples 

We now proceed to summarize the algorithm to solve the optimization prob¬ 
lem. First, let a class of graphs G{a,e) be hxed with the corresponding trans¬ 
lation vectors ei and 62 . Further, the constants Sk£ are constructed by the 
scheme described at the end of the previous section. The main numerical 
step is solution to the uniquely solvable system of linear algebraic equations 
fl3.1ip to construct vertices A and the corresponding double periodic De¬ 
launay graph (AjF). This graph is called optimal in the class G{a,e)- The 
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optimal graph not necessary does correspond to a Voronoi tessellation. In 
this case, one can change a class of graphs by introduction of the new Voronoi 
tessellation for the vertices A. Then, the set A will not necessary be optimal 
in the new class Q'p^. Let (A', i?') be the optimal graph in the class Q'p^. Next, 
if the graph (A', E’') does not correspond to a Voronoi tessellation, it can be 
’’improved” by (A", i?") etc. Therefore, we arrive at the graph chain 

(A, E) (A', E') (A", E”) ^ . (4.1) 

Example 4.1. Consider the hexagonal lattice defined by the fundamental 
translation vectors ei = ^/|(1,0) and 62 = ^)- The area of the cell 

Qq holds unit. Consider N = 3 points (1.075,0.175), (0.919,0.553), (0.444,0.169) 
and the corresponding double periodic Voronoi tessellation shown in Fig|3^. 
Application of the algorithm yields the optimal hexagonal structure FigJSjD. 



Figure 3: a) Three points in the cell Qo are distinguished. Dashed lines show 
the lattice, solid lines the double periodic Delaunay graph, b) The optimal 
graph isomorphic to the graph from a). 


Example 4.2. Consider the hexagonal lattice as in Example 14.II and A = 16 
points with the corresponding double periodic Voronoi tessellation shown in 
Figj2j The considered structure determines a double periodic graph (A, E). 
This graph generates the class of isomorphic graphs G{a,e)- Find the op¬ 
timal graph {A',E) in the class G(a,e)- Construct the Voronoi tessellation 
corresponding to the set A' and the corresponding graph (A', E') which de¬ 
termines the new class The optimal graph in the class is the 

graph {A',E'). Therefore, in this example the graph {A,E) from Fig 12] is 
transformed into the graph (A', E') from FigJU 






















Figure 4: The optimal graph (A', E') from Example I4.21 


Every edge of the Delaunay graph models the interaction caused by the 
diffusion flux. This flux between two disks can be insignihcant if the gap 
between the disk is sufficiently large. In this case, the corresponding edge 
should be deleted from the Delaunay graph. We consider such a case in the 
following example. 

Example 4.3. Consider N = 9 points and the double periodic Voronoi 
tessellation isomorphic to the hexagonal lattice. The perfect hexagonal array 
in FigJS^ presents the optimal graph (A, E). Consider another graph (A, E') 
obtained from (A, E) by deletion of the edges connecting the first layer of 
disks with other layers. In this case, the optimal graph becomes similar to 
hexagonal-layered structure displayed in FiglSb. 



Figure 5: The optimal graphs from Example 14.31 The disks in the layers 1 
and 1’, 2 and 2’, 3 and 3’ are the same in the toroidal topology. 


9 









The above examples illustrate few scenarios of the 2D pattern formations. 
Systematic simulations can help to study properties of the optimal graphs. 
The main feature of the proposed method is investigation of the graph struc¬ 
tures by analytical and numerical methods without treatment of PDE. The 
method of structural approximation recalls a finite element method when a 
continuous problem is approximated by a discrete problem. The structural 
approximation is based on a ’’physical discretization” [2], when edges of the 
graph correspond to the most intensive places of the diffusion flux. Further, 
the principle of minimum energy yields a discrete numerical problem as in a 
hnite element method. The method of structural approximation was justihed 
for the p—Laplacian including linear equations in 0.0.0. 

Remark 4.4. Solution to the optimal energy problem yields solution to the 
optimal packing problem for disks 0 . 0 - The corresponding concentration 
i'{Q) of disks attains the maximal value in the class Q(a,e)- The set of opti¬ 
mal graphs includes graphs corresponding to packing constructed by various 
packing protocols. This scheme gives the set of the optimal concentrations 
depending on protocols, i.e., on the class of graphs. Therefore, in order to 
get the set of all optimal packings, it is sufficient to determine the targets 
of the optimal locations fl4.ip . In the above examples, the graph chain fl4.ip 
terminates. We cannot give an example with an infinite graph chain. 
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